//Use this function to calculate the partial derivatives (Table 2 in Martell et al. 2007)
	//and return estimates of bo and kappa conditional on theta.
	//Arguments: fe=fmsy, ce=MSY
	//tau2 is the standard deviation in recruitment deviations, must add to re to deal with bias correction.

//Added by Chris Grandin and Robyn Forrest May 5 2011
//Modified from TINNS martool.cxx. Original code by Steve Martell

//ONLY BEVERTON HOLT CURRENTLY IMPLEMENTED FOR MANAGEMENT ORIENTED MODEL

void init_model_BH(const dvariable& ce,const dvariable& fe, dvariable& ro, dvariable& kappa,
										const dvariable& m, const dvector& age, const dvar_vector& wa,
										const dvar_vector& fec, const dvar_vector& logsel)
{
	//Use this function to calculate the partial derivatives (Table 2 in Martell et al. 2007)
	//and return estimates of bo and kappa conditional on theta.
	//Arguments: fe=fmsy, ce=MSY

	int i;
	int nage=max(age);
	dvariable phib,phif,phiq,dphif_df=0,dphiq_df=0,dlz_df=0;
	//dvar_vector lx=pow(exp(-m),age-1.); //MOVED into loop to make Chris happy
	//lx(nage)/=(1.-exp(-m));
	dvar_vector va(1,nage);
	dvar_vector lx(1,nage);
	dvar_vector lz(1,nage);

	va=mfexp(logsel);
	dvar_vector za=(m+fe*va);
	dvar_vector sa=1.-exp(-za);
	//dvar_vector fec=elem_prod(wa,ma);

	dvar_vector qa=elem_prod(elem_div(va,za),sa);

	lx(1)=1.0;
	lz(1)=1.0; dlz_df=0;

	for(i=1; i<=nage; i++)
	{
		if(i>1) lx(i)=lx(i-1)*exp(-m);
		if(i>1) lz(i)=lz(i-1)*exp(-za(i-1));
		if(i>1) dlz_df=dlz_df*exp(-za(i-1)) - lz(i-1)*va(i-1)*exp(-za(i-1));
		if(i==nage){ //6/11/2007 added plus group.
					lz(i)/=(1.-mfexp(-za(i)));
					lx(i)/=(1.-mfexp(-m));
					//dlz_df=dlz_df*mfexp(-za(i-1)) - lz(i-1)*va(i-1)*mfexp(-za(i-1))/(1.-mfexp(-za(i)))
					dlz_df=dlz_df/(1.-mfexp(-za(i)))
							-lz(i-1)*mfexp(-za(i-1))*va(i)*mfexp(-za(i))
					/((1.-mfexp(-za(i)))*(1.-mfexp(-za(i))));
				}
		dphif_df=dphif_df+fec(i)*dlz_df;
		dphiq_df=dphiq_df+wa(i)*qa(i)*dlz_df+(lz(i)*wa(i)*va(i)*va(i))/za(i)*(exp(-za(i))-sa(i)/za(i));
	}

	phib=sum(elem_prod(lx,fec)); //unfished eggs
	phif=sum(elem_prod(lz,fec));
	phiq=sum(elem_prod(elem_prod(lz,wa),qa));

	kappa=phib/phif-(fe*phiq*phib/square(phif)*dphif_df)/(phiq+fe*dphiq_df);
	dvariable re=ce/(fe*phiq); //tau2 is the standard deviation in recruitment deviations, must add to re to deal with bias correction??
	ro=re*(kappa-1.)/(kappa-phib/phif);

	//cout<<"Compensation ratio (kappa), Steepness and Ro"<<endl;
	//cout<<kappa<<endl<<kappa/(4+kappa)<<endl<<ro<<endl<<endl;
	//exit(1);
}
